{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "One of the key features of PorePy is the forward-mode automatic (or algorithmic) differentiation, or AD.\n",
    "The AD is used to linearize a system of partial differential equations.\n",
    "\n",
    "This tutorial introduces the basic usage of the AD and shows how it is incorporated into the PorePy models.\n",
    "\n",
    "The time stepping algorithm at the core of PorePy involves numerically solving a discretized system of PDEs $F$ with respect to the vector of variables $x$. \n",
    "Most numerical solution algorithms rely on computing the Jacobian matrix of this system: $J:=\\dfrac{\\partial F}{\\partial x}$. \n",
    "However, programming a formula to compute the Jacobian matrix can be difficult, as it requires analytic differentiation of the equation $F$.\n",
    "\n",
    "Fortunately, the PorePy model developers do not need to perform this task, as the AD framework handles it for them."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AD basics\n",
    "\n",
    "The core of the AD framework is the `AdArray` class. It contains 2 attributes:\n",
    "- `val` - the value of the variable it represents.\n",
    "- `jac` - the Jacobian matrix of the variable it represents.\n",
    "\n",
    "Let's create an independent variable - the vector `x = [1, 2, 3]`.\n",
    "As the variable is independent, its Jacobian is the identity matrix:\n",
    "$\\dfrac{\\partial x}{\\partial x} = I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/firedrake/workspace/porepy_workspace/porepy/src/porepy/numerics/nonlinear/nonlinear_solvers.py:13: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)\n",
      "  from tqdm.autonotebook import trange  # type: ignore\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Values:\n",
      " [1. 2. 3.]\n",
      "Jacobian:\n",
      " [[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sps\n",
    "\n",
    "from porepy.numerics.ad.forward_mode import AdArray\n",
    "\n",
    "val = np.array([1, 2, 3])\n",
    "jac = sps.eye(m=val.size)\n",
    "\n",
    "x = AdArray(val=val, jac=jac)\n",
    "\n",
    "print(\"Values:\\n\", x.val)\n",
    "print(\"Jacobian:\\n\", x.jac.toarray())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use arithmetic operators to build any expression with $x$. For example:\n",
    "$y(x) = 6 x^3 + 3 x - 6$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = 6 * x**3 + 3 * x - 6"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values of the Jacobian $\\dfrac{\\partial y}{\\partial x}$ are available automatically, and they are computed with respect to $x = [1, 2, 3]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 21.   0.   0.]\n",
      " [  0.  75.   0.]\n",
      " [  0.   0. 165.]]\n"
     ]
    }
   ],
   "source": [
    "print(y.jac.toarray())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to defining functions ourselves, we can access a collection of pre-defined standard functions. \n",
    "This collection contains functions such as `sin`, `cos`, `exp` and so on. \n",
    "We refer to the source code for seeing which functions are available, and continue with a simple example of how to access such a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.54030231  0.          0.        ]\n",
      " [ 0.         -0.41614684  0.        ]\n",
      " [ 0.          0.         -0.9899925 ]]\n"
     ]
    }
   ],
   "source": [
    "import porepy.numerics.ad.functions as af\n",
    "\n",
    "y = af.sin(x)\n",
    "print(y.jac.toarray())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The AD framework supports NumPy arrays and SciPy sparse matrices. \n",
    "There are some limitations, but as this will not be covered in this tutorial, the user should check the documentation of the `AdArray` code to become aware of them.\n",
    "\n",
    "Now we are ready to see how the AD is utilized in the PorePy models.\n",
    "\n",
    "# AD in the models"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now consider an example of `SinglePhaseFlow` model, which should be familiar from the [earlier tutorial](./single_phase_flow.ipynb) covering single physics.\n",
    "This time we will not define the geometry manually, but instead use one of the pre-defined geometries provided by a mixin class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import inspect\n",
    "import porepy as pp\n",
    "from porepy.applications.md_grids.model_geometries import (\n",
    "    SquareDomainOrthogonalFractures,\n",
    ")\n",
    "from porepy.models.fluid_mass_balance import SinglePhaseFlow\n",
    "\n",
    "\n",
    "class SinglePhaseFlowWithGeometry(SquareDomainOrthogonalFractures, SinglePhaseFlow):\n",
    "    def meshing_arguments(self) -> dict:\n",
    "        cell_size = self.solid.convert_units(0.25, \"m\")\n",
    "        return {\"cell_size\": cell_size}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Among many other things, PorePy model classes include methods for defining:\n",
    "- the independent variables of the simulation\n",
    "- the system of PDEs describing the physical laws of the simulation\n",
    "\n",
    "All of these are called when the model is initialized in the `prepare_simulation` method.\n",
    "\n",
    "<b>Note:</b>  When running the simulation, you should not call `prepare_simulation` directly. \n",
    "It is called inside `pp.run_time_dependent_model` or `pp.run_stationary_model`. \n",
    "Here, we run it for illustration purposes only, as we do not plan to run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The grid:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHHCAYAAABTMjf2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvR0lEQVR4nO3df1SVVb7H8Q+gnIMpqKmAhjKZ+SN/kHjlopHpUDia5aw1ZdooOmaTP+aqzK00U3QsscYc5iZpmppZjhbXXC0lzbiSmXSdVG6WppU/awR/JRg6kLDvH5NnYkQB5ZwH2O/XWs8i9tn72d+HzSOfnvOcc/yMMUYAAAAW8ne6AAAAAKcQhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEUCfcdddduuuuu5wuo8peffVV+fn56fDhw16fa+bMmfLz86uw38iRIxUZGen1euoSPz8/zZw50+kyvKa2/05c73lW248fV0cQgiMu/cPk5+enbdu2Xfa4MUYRERHy8/PTvffe60CFAAAb1HO6ANjN7XZr1apVuuOOO8q0f/DBB/rmm2/kcrkqtZ/33nvPG+VZacmSJSotLXW6jFrlwoULqlePf05rquHDh+uhhx6q9L8nsAtXhOCoAQMG6K233tLFixfLtK9atUrR0dEKCwur1H4CAwMVGBjojRKtU79+ff5gVJHb7SYIOaywsPCKjwUEBMjtdlfqqWHYhyAERw0dOlSnT5/W5s2bPW3FxcVKT0/XsGHDKr2fytwjtHz5cvn5+WnZsmVl2ufMmSM/Pz9lZGRUOM+7776ruLg43XDDDWrUqJEGDhyozz//vFI1fv755+rXr5+CgoJ000036ZlnnvHalZdt27bp3/7t3+R2u9W2bVu9/PLLlR5b2fshIiMjde+99yorK0s9evRQUFCQunTpoqysLEnS2rVr1aVLF7ndbkVHR2v37t2Vmv/s2bOaNGmSIiIi5HK5dMstt+i5556r1M/qUk3btm1Tz5495Xa7dfPNN+u1116r1NyrV69WdHS0GjVqpODgYHXp0kV//vOfKxxX2XuETpw4odGjRys0NFRut1vdunXTihUrKlXb9R7bvHnz1KtXL914440KCgpSdHS00tPTKzXWGz799FP16dOnzPlw6Ryt6F6ekSNHqmHDhvr66681YMAANWrUSA8//PAV+1flHqF169apc+fOcrvd6ty5s95+++0qHhlqHQM4YPny5UaS+etf/2p69eplhg8f7nls3bp1xt/f33z77bemTZs2ZuDAgRXur0+fPqZPnz4V9rv33ntNSEiIOXr0qDHGmE8//dQEBgaa0aNHVzj2tddeM35+fqZ///7mxRdfNM8995yJjIw0jRs3NocOHbrq2OPHj5vmzZubJk2amJkzZ5o//vGPpl27dqZr165GUoXjq+LTTz81QUFBpnXr1iYlJcXMnj3bhIaGeuaqSGJiomnTpk2F/dq0aWPat29vwsPDzcyZM82f/vQn06pVK9OwYUPz+uuvm9atW5u5c+eauXPnmpCQEHPLLbeYkpKSq+6zsLDQdO3a1dx4443mqaeeMosWLTIjRowwfn5+ZuLEiZWuKTQ01Dz11FNmwYIFpnv37sbPz8989tlnVx373nvvGUnm5z//uUlLSzNpaWlmwoQJ5oEHHqhwXkkmOTn5qn3Onz9vOnbsaOrXr28mT55s/uu//svExcUZSSY1NdWrx2aMMTfddJMZN26cWbBggZk/f77p2bOnkWTWr19f4djK/k5U1jfffGOaNm1qbrzxRjNr1iwzb94806FDB9OtW7dKnQ+JiYnG5XKZtm3bmsTERLNo0SLz2muvXbH/pX9vKtrvpk2bjL+/v+ncubOZP3++mTZtmgkJCTG33XZbtR4/ahaCEBzx0yC0YMEC06hRI3P+/HljjDEPPPCA6du3rzHGVHsQOn78uGnatKm5++67TVFRkbn99ttN69atTX5+/lXHnTt3zjRu3NiMGTOmTHtubq4JCQm5rP1fTZo0yUgy//u//+tpO3HihAkJCan2IDR48GDjdrvNkSNHPG179+41AQEB1R6EJJnt27d72jZt2mQkmaCgoDLzv/zyy0aS2bJly1X3OXv2bHPDDTeYAwcOlGmfMmWKCQgI8ATYimraunWrp+3EiRPG5XKZ3//+91cdO3HiRBMcHGwuXrx41X7lqUwQSk1NNZLM66+/7mkrLi42sbGxpmHDhqagoOCq46/n2IwxnvPrp3N37tzZ9OvXr8Kx1R2Efve73xk/Pz+ze/duT9vp06dN06ZNKx2EJJkpU6ZUar7KBqGoqCgTHh5uzp4962m7FJAJQnUXT43BcQ8++KAuXLig9evX69y5c1q/fn2VnharirCwMKWlpWnz5s2Ki4tTTk6Oli1bpuDg4KuO27x5s86ePauhQ4fq1KlTni0gIEAxMTHasmXLVcdnZGTo3//939WzZ09PW/Pmza96Of9alJSUaNOmTRo8eLBat27tae/YsaMSEhKqdS5J6tSpk2JjYz3fx8TESJL69etXZv5L7QcPHrzq/t566y3FxcWpSZMmZX7O8fHxKikp0datWytVU1xcnOf75s2bq3379hXO3bhxYxUWFpZ5mrY6ZWRkKCwsTEOHDvW01a9fX//xH/+h77//Xh988EGF+7jWY5OkoKAgz39/9913ys/PV1xcnHbt2lXFI7l+GzduVGxsrKKiojxtTZs2rfL5MHbs2Gqr6fjx48rJyVFiYqJCQkI87Xfffbc6depUbfOg5uHuPjiuefPmio+P16pVq3T+/HmVlJToV7/6ldfme+ihh/T6669rw4YNevTRR/Xzn/+8wjFffvmlpH/8gS9PRUHqyJEjnjDwU+3bt69w7gsXLig/P79M25VuIj958qQuXLigdu3alTtXZe6Dqoqfhh1Jnj8gERER5bZ/9913V93fl19+qU8//VTNmzcv9/ETJ05UuSZJatKkSYVzjxs3Tm+++aZ+8YtfqFWrVrrnnnv04IMPqn///hXOWRlHjhxRu3bt5O9f9v8/O3bs6Hm8Itd6bJK0fv16PfPMM8rJyVFRUZGn3YkbiI8cOVImQF9yyy23VHof9erV00033VStNUm64rnjRGCEbxCEUCMMGzZMY8aMUW5urn7xi1+ocePGXpvr9OnT+uSTTyRJe/fuVWlp6WV/nP7VpRt1V65cWW4I8eYrhtasWaNRo0aVaTPGeG2+qggICKhSe0V1l5aW6u6779YTTzxR7uO33nrrNddU0dwtWrRQTk6ONm3apHfffVfvvvuuli9frhEjRlT6hmZvu9Zj+/DDD3Xffffpzjvv1EsvvaTw8HDVr19fy5cv16pVq7xRqte5XK4Kz1ugMghCqBF++ctf6re//a0+/vhjrVmzxqtzjR8/XufOnVNKSoqmTp2q1NRUJSUlXXVM27ZtJf3jj2V8fHyV52zTpo3nqtJP7d+/v8KxCQkJlX66pnnz5goKCrrmuZzWtm1bff/999f0M64OgYGBGjRokAYNGqTS0lKNGzdOL7/8sqZPn16lqxXladOmjT799NPLgvcXX3zhedxb/vu//1tut1ubNm0q89YIy5cv99qcV9OmTRt99dVXl7WX1+Yrl37+tfXcwbUjTqNGaNiwoRYuXKiZM2dq0KBBXpsnPT1da9as0dy5czVlyhQ99NBDevrpp3XgwIGrjktISFBwcLDmzJmjH3744bLHT548edXxAwYM0Mcff6wdO3aUGfPGG29UWHN4eLji4+PLbFcSEBCghIQErVu3TkePHvW079u3T5s2bapwLqc9+OCDys7OLrfWs2fPXvZ+U9Xp9OnTZb739/dX165dJanMU0nXasCAAcrNzS0T9C9evKgXX3xRDRs2VJ8+fa57jisJCAiQn5+fSkpKPG2HDx/WunXrvDbn1SQkJCg7O1s5OTmetjNnzlTqfPCW8PBwRUVFacWKFWWeit68ebP27t3rWF3wPq4IocZITEz06v5PnDihsWPHqm/fvpowYYIkacGCBdqyZYtGjhypbdu2XfFSe3BwsBYuXKjhw4ere/fueuihh9S8eXMdPXpUGzZsUO/evbVgwYIrzv3EE09o5cqV6t+/vyZOnKgbbrhBixcv9lwlqE6zZs3Sxo0bFRcXp3Hjxnn+2N52223VPld1e/zxx/XOO+/o3nvv1ciRIxUdHa3CwkLt2bNH6enpOnz4sJo1a+aVuR955BGdOXNG/fr100033aQjR47oxRdfVFRUlOc+nuvx6KOP6uWXX9bIkSO1c+dORUZGKj09XR999JFSU1PVqFGjajiK8g0cOFDz589X//79NWzYMJ04cUJpaWm65ZZbHPmdeOKJJ/T666/r7rvv1u9+9zvdcMMNeuWVV9S6dWudOXPGsTc+TElJ0cCBA3XHHXfoN7/5jc6cOeM5d77//ntHaoL3cUUI1hg7dqyKioo8b9omSTfeeKMWL16s7OxszZs376rjhw0bpszMTLVq1Up//OMfNXHiRK1evVpRUVGX3cPzr8LDw7VlyxZ17dpVc+fOVWpqqkaMGKGJEydW2/Fd0rVrV23atEnNmzfXjBkztGzZMs2aNUu//OUvq32u6tagQQN98MEHevzxx5WVlaWJEydq7ty5+vLLLzVr1qwyr+apbr/+9a/ldrv10ksvady4cVqxYoWGDBmid999t1ruRQkKClJWVpYefvhhrVixQr///e915swZLV++3Cu/Bz/Vr18/LV26VLm5uZo0aZL+8pe/6LnnnnPsdyIiIkJbtmxRx44dNWfOHKWmpioxMVG/+c1vJP3jnbqd0L9/f7311lsqKSnR1KlTtXbtWi1fvlw9evRwpB74hp+pKXddAgCsNmnSJL388sv6/vvvr3hjOFDduCIEAPC5CxculPn+9OnTWrlype644w5CEHyKe4QAAD4XGxuru+66Sx07dlReXp6WLl2qgoICTZ8+3enSYBmCEADA5wYMGKD09HQtXrxYfn5+6t69u5YuXao777zT6dJgGUefGtu6dasGDRqkli1bys/Pr1Iv5czKylL37t09n0r96quver1OAED1mjNnjg4cOKDz58+rsLBQH374oWPvHwW7ORqECgsL1a1bN6WlpVWq/6FDhzRw4ED17dtXOTk5mjRpkh555JFa8f4oAACg5qkxrxrz8/PT22+/rcGDB1+xz5NPPqkNGzbos88+87Q99NBDOnv2rDZu3OiDKgEAQF1Sq+4Rys7OvuzSaUJCgiZNmnTFMUVFRWXeFba0tFRnzpzRjTfe6NibdgEAgKoxxujcuXNq2bJltX7OXK0KQrm5uQoNDS3TFhoaqoKCAl24cEFBQUGXjUlJSdGsWbN8VSIAAPCiY8eO6aabbqq2/dWqIHQtpk6dWuYDNfPz89W6dWv9WVKUY1XBV96T9KykxZLaO1wLvO81SUslTZN0j8O1wPs4v+2SI2miVO0fR1OrglBYWJjy8vLKtOXl5Sk4OLjcq0GS5HK5ynza8iVRkniRZt137Mev0ZK6O1kIfGLrj187ivPbBpzfdqru21pq1TtLx8bGKjMzs0zb5s2bFRsb61BFAACgNnM0CH3//ffKyclRTk6OpH+8PD4nJ0dHjx6V9I+ntUaMGOHp/9hjj+ngwYN64okn9MUXX+ill17Sm2++qcmTJztRPgAAqOUcDUKffPKJbr/9dt1+++2SpKSkJN1+++2aMWOGJOn48eOeUCRJP/vZz7RhwwZt3rxZ3bp10wsvvKBXXnlFCQkJjtQPAABqN0fvEbrrrrt0tbcxKu9do++66y7t3r3bi1UBAABb1Kp7hAAAAKoTQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWo4HobS0NEVGRsrtdismJkY7duy4av/U1FS1b99eQUFBioiI0OTJk/X3v//dR9UCAIC6xNEgtGbNGiUlJSk5OVm7du1St27dlJCQoBMnTpTbf9WqVZoyZYqSk5O1b98+LV26VGvWrNFTTz3l48oBAEBd4GgQmj9/vsaMGaNRo0apU6dOWrRokRo0aKBly5aV23/79u3q3bu3hg0bpsjISN1zzz0aOnRohVeRAAAAyuNYECouLtbOnTsVHx//z2L8/RUfH6/s7Oxyx/Tq1Us7d+70BJ+DBw8qIyNDAwYMuOI8RUVFKigoKLMBAABIUj2nJj516pRKSkoUGhpapj00NFRffPFFuWOGDRumU6dO6Y477pAxRhcvXtRjjz121afGUlJSNGvWrGqtHQAA1A2O3yxdFVlZWZozZ45eeukl7dq1S2vXrtWGDRs0e/bsK46ZOnWq8vPzPduxY8d8WDEAAKjJHLsi1KxZMwUEBCgvL69Me15ensLCwsodM336dA0fPlyPPPKIJKlLly4qLCzUo48+qmnTpsnf//Jc53K55HK5qv8AAABArefYFaHAwEBFR0crMzPT01ZaWqrMzEzFxsaWO+b8+fOXhZ2AgABJkjHGe8UCAIA6ybErQpKUlJSkxMRE9ejRQz179lRqaqoKCws1atQoSdKIESPUqlUrpaSkSJIGDRqk+fPn6/bbb1dMTIy++uorTZ8+XYMGDfIEIgAAgMpyNAgNGTJEJ0+e1IwZM5Sbm6uoqCht3LjRcwP10aNHy1wBevrpp+Xn56enn35a3377rZo3b65Bgwbp2WefdeoQAABALeZnLHtOqaCgQCEhIfpA0p1OFwOve0PSryXtlNTd4Vrgfc9Imi7pdUkPO1wLvI/z2y5bJfWRlJ+fr+Dg4Grbb6161RgAAEB1IggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUc/dBVJ+2X1NDpIuB1h378us/RKuArf/vx6yFJu5wsBD7B+W2X/V7ar7UfugoAAGqf6v7QVWuvCE2T1NHpIuB1H0laKGm2pJ85XAu8b52kdEljJfV2thT4AOe3XfZJetYL+7U2CN0j6U6ni4BPLJQ0QFJ3pwuB1x3SP4JQb0kPO1wLfIPz2x5b5Z0gxM3SAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtRwPQmlpaYqMjJTb7VZMTIx27Nhx1f5nz57V+PHjFR4eLpfLpVtvvVUZGRk+qhYAANQl9ZycfM2aNUpKStKiRYsUExOj1NRUJSQkaP/+/WrRosVl/YuLi3X33XerRYsWSk9PV6tWrXTkyBE1btzY98UDAIBaz9EgNH/+fI0ZM0ajRo2SJC1atEgbNmzQsmXLNGXKlMv6L1u2TGfOnNH27dtVv359SVJkZKQvSwYAAHWIY0+NFRcXa+fOnYqPj/9nMf7+io+PV3Z2drlj3nnnHcXGxmr8+PEKDQ1V586dNWfOHJWUlFxxnqKiIhUUFJTZAAAAJAeD0KlTp1RSUqLQ0NAy7aGhocrNzS13zMGDB5Wenq6SkhJlZGRo+vTpeuGFF/TMM89ccZ6UlBSFhIR4toiIiGo9DgAAUHs5frN0VZSWlqpFixZavHixoqOjNWTIEE2bNk2LFi264pipU6cqPz/fsx07dsyHFQMAgJrMsXuEmjVrpoCAAOXl5ZVpz8vLU1hYWLljwsPDVb9+fQUEBHjaOnbsqNzcXBUXFyswMPCyMS6XSy6Xq3qLBwAAdYJjV4QCAwMVHR2tzMxMT1tpaakyMzMVGxtb7pjevXvrq6++UmlpqaftwIEDCg8PLzcEAQAAXI2jT40lJSVpyZIlWrFihfbt26exY8eqsLDQ8yqyESNGaOrUqZ7+Y8eO1ZkzZzRx4kQdOHBAGzZs0Jw5czR+/HinDgEAANRijr58fsiQITp58qRmzJih3NxcRUVFaePGjZ4bqI8ePSp//39mtYiICG3atEmTJ09W165d1apVK02cOFFPPvmkU4cAAABqMUeDkCRNmDBBEyZMKPexrKysy9piY2P18ccfe7kqAABgg1r1qjEAAIDqRBACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrVTkIJSYmauvWrd6oBQAAwKeqHITy8/MVHx+vdu3aac6cOfr222+9URcAAIDX1avqgHXr1unkyZNauXKlVqxYoeTkZMXHx2v06NG6//77Vb9+fW/UWe32S2rodBHwukM/ft3naBXwlb/9+PWQpF1OFgKf4Py2y34v7dfPGGOuZwe7du3S8uXL9corr6hhw4b69a9/rXHjxqldu3bVVWO1KigoUEhIiNNlAACAa5Cfn6/g4OBq21+Vrwj91PHjx7V582Zt3rxZAQEBGjBggPbs2aNOnTrp+eef1+TJk6urzmo3TVJHp4uA130kaaGk2ZJ+5nAt8L51ktIljZXU29lS4AOc33bZJ+lZb+zYVFFxcbFJT083AwcONPXr1zfR0dFm4cKFJj8/39Nn7dq1pnHjxlXdtU/k5+cbSeYDyRi2Or+9LhlJZmcNqIXN+9vsH9f79RpQC5v3N85vu7YPflzvn+aN6lDlK0Lh4eEqLS3V0KFDtWPHDkVFRV3Wp2/fvmrcuPH1JTQAAAAvq3IQ+tOf/qQHHnhAbrf7in0aN26sQ4cOXVdhAAAA3lblIDR8+HBv1AEAAOBzvLM0AACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsVSOCUFpamiIjI+V2uxUTE6MdO3ZUatzq1avl5+enwYMHe7dAAABQJzkehNasWaOkpCQlJydr165d6tatmxISEnTixImrjjt8+LD+8z//U3FxcT6qFAAA1DWOB6H58+drzJgxGjVqlDp16qRFixapQYMGWrZs2RXHlJSU6OGHH9asWbN08803+7BaAABQlzgahIqLi7Vz507Fx8d72vz9/RUfH6/s7OwrjvvDH/6gFi1aaPTo0RXOUVRUpIKCgjIbAACA5HAQOnXqlEpKShQaGlqmPTQ0VLm5ueWO2bZtm5YuXaolS5ZUao6UlBSFhIR4toiIiOuuGwAA1A2OPzVWFefOndPw4cO1ZMkSNWvWrFJjpk6dqvz8fM927NgxL1cJAABqi3pOTt6sWTMFBAQoLy+vTHteXp7CwsIu6//111/r8OHDGjRokKettLRUklSvXj3t379fbdu2LTPG5XLJ5XJ5oXoAAFDbOXpFKDAwUNHR0crMzPS0lZaWKjMzU7GxsZf179Chg/bs2aOcnBzPdt9996lv377KycnhaS8AAFAljl4RkqSkpCQlJiaqR48e6tmzp1JTU1VYWKhRo0ZJkkaMGKFWrVopJSVFbrdbnTt3LjO+cePGknRZOwAAQEUcD0JDhgzRyZMnNWPGDOXm5ioqKkobN2703EB99OhR+fvXqluZAABALeF4EJKkCRMmaMKECeU+lpWVddWxr776avUXBAAArMClFgAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLXqOV2AU/ZLauh0EfC6Qz9+3edoFfCVv/349ZCkXU4WAp/g/LbLfi/t188YY7y07xqpoKBAISEhTpcBAACuQX5+voKDg6ttf9ZeEZomqaPTRcDrPpK0UNJsST9zuBZ43zpJ6ZLGSurtbCnwAc5vu+yT9KwX9mttELpH0p1OFwGfWChpgKTuThcCrzukfwSh3pIedrgW+Abntz22yjtBiJulAQCAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALBWjQhCaWlpioyMlNvtVkxMjHbs2HHFvkuWLFFcXJyaNGmiJk2aKD4+/qr9AQAArsTxILRmzRolJSUpOTlZu3btUrdu3ZSQkKATJ06U2z8rK0tDhw7Vli1blJ2drYiICN1zzz369ttvfVw5AACo7RwPQvPnz9eYMWM0atQoderUSYsWLVKDBg20bNmycvu/8cYbGjdunKKiotShQwe98sorKi0tVWZmpo8rBwAAtZ2jQai4uFg7d+5UfHy8p83f31/x8fHKzs6u1D7Onz+vH374QU2bNi338aKiIhUUFJTZAAAAJIeD0KlTp1RSUqLQ0NAy7aGhocrNza3UPp588km1bNmyTJj6qZSUFIWEhHi2iIiI664bAADUDY4/NXY95s6dq9WrV+vtt9+W2+0ut8/UqVOVn5/v2Y4dO+bjKgEAQE1Vz8nJmzVrpoCAAOXl5ZVpz8vLU1hY2FXHzps3T3PnztX777+vrl27XrGfy+WSy+WqlnoBAEDd4ugVocDAQEVHR5e50fnSjc+xsbFXHPf8889r9uzZ2rhxo3r06OGLUgEAQB3k6BUhSUpKSlJiYqJ69Oihnj17KjU1VYWFhRo1apQkacSIEWrVqpVSUlIkSc8995xmzJihVatWKTIy0nMvUcOGDdWwYUPHjgMAANQ+jgehIUOG6OTJk5oxY4Zyc3MVFRWljRs3em6gPnr0qPz9/3nhauHChSouLtavfvWrMvtJTk7WzJkzfVk6AACo5RwPQpI0YcIETZgwodzHsrKyynx/+PBh7xcEAACsUKtfNQYAAHA9CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArEUQAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBaBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEAAAsBZBCAAAWIsgBAAArFUjglBaWpoiIyPldrsVExOjHTt2XLX/W2+9pQ4dOsjtdqtLly7KyMjwUaUAAKAucTwIrVmzRklJSUpOTtauXbvUrVs3JSQk6MSJE+X23759u4YOHarRo0dr9+7dGjx4sAYPHqzPPvvMx5UDAIDazvEgNH/+fI0ZM0ajRo1Sp06dtGjRIjVo0EDLli0rt/+f//xn9e/fX48//rg6duyo2bNnq3v37lqwYIGPKwcAALWdo0GouLhYO3fuVHx8vKfN399f8fHxys7OLndMdnZ2mf6SlJCQcMX+AAAAV1LPyclPnTqlkpIShYaGlmkPDQ3VF198Ue6Y3Nzccvvn5uaW27+oqEhFRUWe7/Pz8yVJOddRN2qPfT9+3SnpeycLgU8c/vHrPklbHawDvsH5bZecH78aY6p1v44GIV9ISUnRrFmzLmuf6EAtcM6jThcAn3r2xw124Py2y+nTpxUSElJt+3M0CDVr1kwBAQHKy8sr056Xl6ewsLByx4SFhVWp/9SpU5WUlOT5/uzZs2rTpo2OHj1arT9IVF1BQYEiIiJ07NgxBQcHO12O9ViPmoO1qDlYi5ojPz9frVu3VtOmTat1v44GocDAQEVHRyszM1ODBw+WJJWWliozM1MTJkwod0xsbKwyMzM1adIkT9vmzZsVGxtbbn+XyyWXy3VZe0hICL/UNURwcDBrUYOwHjUHa1FzsBY1h79/9d7e7PhTY0lJSUpMTFSPHj3Us2dPpaamqrCwUKNGjZIkjRgxQq1atVJKSookaeLEierTp49eeOEFDRw4UKtXr9Ynn3yixYsXO3kYAACgFnI8CA0ZMkQnT57UjBkzlJubq6ioKG3cuNFzQ/TRo0fLpL9evXpp1apVevrpp/XUU0+pXbt2WrdunTp37uzUIQAAgFrK8SAkSRMmTLjiU2FZWVmXtT3wwAN64IEHrmkul8ul5OTkcp8ug2+xFjUL61FzsBY1B2tRc3hrLfxMdb8ODQAAoJZw/J2lAQAAnEIQAgAA1iIIAQAAaxGEAACAtepkEEpLS1NkZKTcbrdiYmK0Y8eOq/Z/66231KFDB7ndbnXp0kUZGRk+qrTuq8paLFmyRHFxcWrSpImaNGmi+Pj4CtcOVVPVc+OS1atXy8/Pz/PGp7h+VV2Ls2fPavz48QoPD5fL5dKtt97Kv1XVpKprkZqaqvbt2ysoKEgRERGaPHmy/v73v/uo2rpr69atGjRokFq2bCk/Pz+tW7euwjFZWVnq3r27XC6XbrnlFr366qtVn9jUMatXrzaBgYFm2bJl5vPPPzdjxowxjRs3Nnl5eeX2/+ijj0xAQIB5/vnnzd69e83TTz9t6tevb/bs2ePjyuueqq7FsGHDTFpamtm9e7fZt2+fGTlypAkJCTHffPONjyuvm6q6HpccOnTItGrVysTFxZn777/fN8XWcVVdi6KiItOjRw8zYMAAs23bNnPo0CGTlZVlcnJyfFx53VPVtXjjjTeMy+Uyb7zxhjl06JDZtGmTCQ8PN5MnT/Zx5XVPRkaGmTZtmlm7dq2RZN5+++2r9j948KBp0KCBSUpKMnv37jUvvviiCQgIMBs3bqzSvHUuCPXs2dOMHz/e831JSYlp2bKlSUlJKbf/gw8+aAYOHFimLSYmxvz2t7/1ap02qOpa/KuLFy+aRo0amRUrVnirRKtcy3pcvHjR9OrVy7zyyismMTGRIFRNqroWCxcuNDfffLMpLi72VYnWqOpajB8/3vTr169MW1JSkundu7dX67RNZYLQE088YW677bYybUOGDDEJCQlVmqtOPTVWXFysnTt3Kj4+3tPm7++v+Ph4ZWdnlzsmOzu7TH9JSkhIuGJ/VM61rMW/On/+vH744Ydq/4A9G13revzhD39QixYtNHr0aF+UaYVrWYt33nlHsbGxGj9+vEJDQ9W5c2fNmTNHJSUlviq7TrqWtejVq5d27tzpefrs4MGDysjI0IABA3xSM/6puv5+14h3lq4up06dUklJiefjOS4JDQ3VF198Ue6Y3Nzccvvn5uZ6rU4bXMta/Ksnn3xSLVu2vOwXHVV3Leuxbds2LV26VDk5OT6o0B7XshYHDx7U//zP/+jhhx9WRkaGvvrqK40bN04//PCDkpOTfVF2nXQtazFs2DCdOnVKd9xxh4wxunjxoh577DE99dRTvigZP3Glv98FBQW6cOGCgoKCKrWfOnVFCHXH3LlztXr1ar399ttyu91Ol2Odc+fOafjw4VqyZImaNWvmdDnWKy0tVYsWLbR48WJFR0dryJAhmjZtmhYtWuR0adbJysrSnDlz9NJLL2nXrl1au3atNmzYoNmzZztdGq5Rnboi1KxZMwUEBCgvL69Me15ensLCwsodExYWVqX+qJxrWYtL5s2bp7lz5+r9999X165dvVmmNaq6Hl9//bUOHz6sQYMGedpKS0slSfXq1dP+/fvVtm1b7xZdR13LuREeHq769esrICDA09axY0fl5uaquLhYgYGBXq25rrqWtZg+fbqGDx+uRx55RJLUpUsXFRYW6tFHH9W0adPKfEg4vOtKf7+Dg4MrfTVIqmNXhAIDAxUdHa3MzExPW2lpqTIzMxUbG1vumNjY2DL9JWnz5s1X7I/KuZa1kKTnn39es2fP1saNG9WjRw9flGqFqq5Hhw4dtGfPHuXk5Hi2++67T3379lVOTo4iIiJ8WX6dci3nRu/evfXVV195wqgkHThwQOHh4YSg63Ata3H+/PnLws6lgGr46E6fqra/31W7j7vmW716tXG5XObVV181e/fuNY8++qhp3Lixyc3NNcYYM3z4cDNlyhRP/48++sjUq1fPzJs3z+zbt88kJyfz8vlqUtW1mDt3rgkMDDTp6enm+PHjnu3cuXNOHUKdUtX1+Fe8aqz6VHUtjh49aho1amQmTJhg9u/fb9avX29atGhhnnnmGacOoc6o6lokJyebRo0amb/85S/m4MGD5r333jNt27Y1Dz74oFOHUGecO3fO7N692+zevdtIMvPnzze7d+82R44cMcYYM2XKFDN8+HBP/0svn3/88cfNvn37TFpaGi+fv+TFF180rVu3NoGBgaZnz57m448/9jzWp08fk5iYWKb/m2++aW699VYTGBhobrvtNrNhwwYfV1x3VWUt2rRpYyRdtiUnJ/u+8DqqqufGTxGEqldV12L79u0mJibGuFwuc/PNN5tnn33WXLx40cdV101VWYsffvjBzJw507Rt29a43W4TERFhxo0bZ7777jvfF17HbNmypdy/AZd+/omJiaZPnz6XjYmKijKBgYHm5ptvNsuXL6/yvH7GcC0PAADYqU7dIwQAAFAVBCEAAGAtghAAALAWQQgAAFiLIAQAAKxFEAIAANYiCAEAAGsRhAAAgLUIQgAAwFoEIQAAYC2CEIBa7+TJkwoLC9OcOXM8bdu3b1dgYOBln04NAD/FZ40BqBMyMjI0ePBgbd++Xe3bt1dUVJTuv/9+zZ8/3+nSANRgBCEAdcb48eP1/vvvq0ePHtqzZ4/++te/yuVyOV0WgBqMIASgzrhw4YI6d+6sY8eOaefOnerSpYvTJQGo4bhHCECd8fXXX+tvf/ubSktLdfjwYafLAVALcEUIQJ1QXFysnj17KioqSu3bt1dqaqr27NmjFi1aOF0agBqMIASgTnj88ceVnp6u//u//1PDhg3Vp08fhYSEaP369U6XBqAG46kxALVeVlaWUlNTtXLlSgUHB8vf318rV67Uhx9+qIULFzpdHoAajCtCAADAWlwRAgAA1iIIAQAAaxGEAACAtQhCAADAWgQhAABgLYIQAACwFkEIAABYiyAEAACsRRACAADWIggBAABrEYQAAIC1CEIAAMBa/w/+r52ltDro3AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "One fracture inside:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAGiCAYAAACGUJO6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAc7klEQVR4nO3df2zV9b3H8Vd74JwCpQXX9afFBpw/UZhl9FbkMpdOElwdMcZOvLTjKk7tFqWZQOVHVZQiKummVQITNbm6ogSNGU2ddhKjdCErNHETMRWwDHcKzHEKB23h9HP/mBxXW6Df2p7yts9H8k3s1++Pdz+pfXpOT0/jnHNOAAAYEz/YAwAA0BcEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGCS54C98847KiwsVGZmpuLi4vTaa6+d9ZytW7fqqquuUiAQ0IUXXqjnn3++D6MCAPAVzwELh8OaNGmSqqure3X83r17df311+vaa69VU1OT7r33Xt1+++164403PA8LAMApcd/kzXzj4uL06quvavbs2ac9ZtGiRdqyZYv++te/Rvf97Gc/05EjR1RXV9fXWwMAhrhhA32DhoYGFRQUdNk3c+ZM3Xvvvac9p729Xe3t7dGPOzs79dlnn+k73/mO4uLiBmpUAMAAcM7p6NGjyszMVHx8/730YsADFgwGlZaW1mVfWlqa2tra9Pnnn2vEiBHdzqmsrNSDDz440KMBAGJo//79Ov/88/vtegMesL4oLy9XWVlZ9ONQKKRx48Z1O+6B++/XHf/7v7EcDegX2/78Z930P/9z1uM2/d//6er/+q8YTAT0n3UbNuiBlSu77R89enS/3mfAA5aenq7W1tYu+1pbW5WUlNTjoy9JCgQCCgQCZ7yuz+dTeUWF/H5/v80KxMrsnBydv3ixDhw4oJ5+DB0XF6fzzz9fs3/2M/l8vkGYEOi78ooKrXj0UUUikS77+/tHQAP+e2D5+fmqr6/vsu/NN99Ufn7+N7puWVkZ8YJZPp9Pv/nNbyR1/4/61MdVVVXECyb5/f4uz6INFM8BO3bsmJqamtTU1CTp3y+Tb2pqUktLi6R/P/1XXFwcPf7OO+/Unj17tHDhQn344Yd6+umn9fLLL2vBggV9Gtjn8+m+++7T6tWr+3Q+cK648cYbtWnTJmVmZnbZf/7552vTpk268cYbB2ky4JtbvXq17rvvvn590UY3zqO3337bSeq2lZSUOOecKykpcTNmzOh2zuTJk53f73fjx493zz33nKd7hkIhJ8k9cP/9rr293evIwDnt1Ne3JFdbW+tOnjw52CMB/eaz/fujX9+hUKhfr/2Nfg8sVtra2pScnKxPm5uVMWHCYI8D9KtwOKzExERJ/36GY9SoUYM8EdB/wgcPKvHLV6KHQiElJSX127V5L0QAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJjUp4BVV1crJydHCQkJysvL0/bt2894fFVVlS6++GKNGDFC2dnZWrBggb744os+DQwAgNSHgG3cuFFlZWWqqKjQjh07NGnSJM2cOVMHDx7s8fiXXnpJixcvVkVFhXbt2qVnn31WGzdu1P333/+NhwcADF2eA7ZmzRrNnz9f8+bN02WXXaa1a9dq5MiR2rBhQ4/Hb9u2TdOmTdOcOXOUk5Oj6667TrfccstZH7UBAHAmngLW0dGhxsZGFRQUfHWB+HgVFBSooaGhx3OuvvpqNTY2RoO1Z88e1dbWatasWae9T3t7u9ra2rpsAAD8p2FeDj58+LAikYjS0tK67E9LS9OHH37Y4zlz5szR4cOHdc0118g5p5MnT+rOO+8841OIlZWVevDBB72MBgAYYgb8VYhbt27VypUr9fTTT2vHjh3avHmztmzZohUrVpz2nPLycoVCoei2f//+gR4TAGCMp0dgKSkp8vl8am1t7bK/tbVV6enpPZ6zbNkyzZ07V7fffrsk6YorrlA4HNYdd9yhJUuWKD6+e0MDgYACgYCX0QAAQ4ynR2B+v1+5ubmqr6+P7uvs7FR9fb3y8/N7POf48ePdIuXz+SRJzjmv8wIAIMnjIzBJKisrU0lJiaZMmaKpU6eqqqpK4XBY8+bNkyQVFxcrKytLlZWVkqTCwkKtWbNG3//+95WXl6fm5mYtW7ZMhYWF0ZABAOCV54AVFRXp0KFDWr58uYLBoCZPnqy6urroCztaWlq6POJaunSp4uLitHTpUh04cEDf/e53VVhYqEceeaT/PgsAwJAT5ww8j9fW1qbk5GR92tysjAkTBnscoF+Fw2ElJiZKko4dO6ZRo0YN8kRA/wkfPKjELx/ghEIhJSUl9du1eS9EAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCY1KeAVVdXKycnRwkJCcrLy9P27dvPePyRI0dUWlqqjIwMBQIBXXTRRaqtre3TwAAASNIwryds3LhRZWVlWrt2rfLy8lRVVaWZM2dq9+7dSk1N7XZ8R0eHfvzjHys1NVWbNm1SVlaWPvnkE40ZM6Y/5gcADFGeA7ZmzRrNnz9f8+bNkyStXbtWW7Zs0YYNG7R48eJux2/YsEGfffaZtm3bpuHDh0uScnJyvtnUAIAhz9NTiB0dHWpsbFRBQcFXF4iPV0FBgRoaGno85/XXX1d+fr5KS0uVlpamiRMnauXKlYpEIqe9T3t7u9ra2rpsAAD8J08BO3z4sCKRiNLS0rrsT0tLUzAY7PGcPXv2aNOmTYpEIqqtrdWyZcv0xBNP6OGHHz7tfSorK5WcnBzdsrOzvYwJABgCBvxViJ2dnUpNTdW6deuUm5uroqIiLVmyRGvXrj3tOeXl5QqFQtFt//79Az0mAMAYTz8DS0lJkc/nU2tra5f9ra2tSk9P7/GcjIwMDR8+XD6fL7rv0ksvVTAYVEdHh/x+f7dzAoGAAoGAl9EAAEOMp0dgfr9fubm5qq+vj+7r7OxUfX298vPzezxn2rRpam5uVmdnZ3TfRx99pIyMjB7jBQBAb3h+CrGsrEzr16/XCy+8oF27dumuu+5SOByOviqxuLhY5eXl0ePvuusuffbZZ7rnnnv00UcfacuWLVq5cqVKS0v777MAAAw5nl9GX1RUpEOHDmn58uUKBoOaPHmy6urqoi/saGlpUXz8V13Mzs7WG2+8oQULFujKK69UVlaW7rnnHi1atKj/PgsAwJAT55xzgz3E2bS1tSk5OVmfNjcrY8KEwR4H6FfhcFiJiYmSpGPHjmnUqFGDPBHQf8IHDyrxywc4oVBISUlJ/XZt3gsRAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABgEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAm9Slg1dXVysnJUUJCgvLy8rR9+/ZenVdTU6O4uDjNnj27L7cFACDKc8A2btyosrIyVVRUaMeOHZo0aZJmzpypgwcPnvG8ffv26de//rWmT5/e52EBADjFc8DWrFmj+fPna968ebrsssu0du1ajRw5Uhs2bDjtOZFIRLfeeqsefPBBjR8//qz3aG9vV1tbW5cNAID/5ClgHR0damxsVEFBwVcXiI9XQUGBGhoaTnveQw89pNTUVN122229uk9lZaWSk5OjW3Z2tpcxAQBDgKeAHT58WJFIRGlpaV32p6WlKRgM9njOu+++q2effVbr16/v9X3Ky8sVCoWi2/79+72MCQAYAoYN5MWPHj2quXPnav369UpJSen1eYFAQIFAYAAnAwBY5ylgKSkp8vl8am1t7bK/tbVV6enp3Y7/+OOPtW/fPhUWFkb3dXZ2/vvGw4Zp9+7dmjBhQl/mBgAMcZ6eQvT7/crNzVV9fX10X2dnp+rr65Wfn9/t+EsuuUTvv/++mpqaotsNN9yga6+9Vk1NTfxsCwDQZ56fQiwrK1NJSYmmTJmiqVOnqqqqSuFwWPPmzZMkFRcXKysrS5WVlUpISNDEiRO7nD9mzBhJ6rYfAAAvPAesqKhIhw4d0vLlyxUMBjV58mTV1dVFX9jR0tKi+Hje4AMAMLDinHNusIc4m7a2NiUnJ+vT5mZl8DMzfMuEw2ElJiZKko4dO6ZRo0YN8kRA/wkfPKjELx/ghEIhJSUl9du1eagEADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACT+hSw6upq5eTkKCEhQXl5edq+fftpj12/fr2mT5+usWPHauzYsSooKDjj8QAA9IbngG3cuFFlZWWqqKjQjh07NGnSJM2cOVMHDx7s8fitW7fqlltu0dtvv62GhgZlZ2fruuuu04EDB77x8ACAoSvOOee8nJCXl6cf/OAHeuqppyRJnZ2dys7O1q9+9SstXrz4rOdHIhGNHTtWTz31lIqLi3s8pr29Xe3t7dGP29ralJ2drU+bm5UxYYKXcYFzXjgcVmJioiTp2LFjGjVq1CBPBPSf8MGDSkxLkySFQiElJSX127U9PQLr6OhQY2OjCgoKvrpAfLwKCgrU0NDQq2scP35cJ06c0HnnnXfaYyorK5WcnBzdsrOzvYwJABgCPAXs8OHDikQiSvuypqekpaUpGAz26hqLFi1SZmZmlwh+XXl5uUKhUHTbv3+/lzEBAEPAsFjebNWqVaqpqdHWrVuVkJBw2uMCgYACgUAMJwMAWOMpYCkpKfL5fGptbe2yv7W1Venp6Wc89/HHH9eqVav01ltv6corr/Q+KQAA/8HTU4h+v1+5ubmqr6+P7uvs7FR9fb3y8/NPe97q1au1YsUK1dXVacqUKX2fFgCAL3l+CrGsrEwlJSWaMmWKpk6dqqqqKoXDYc2bN0+SVFxcrKysLFVWVkqSHn30US1fvlwvvfSScnJyoj8rS0xMjL7yCgAArzwHrKioSIcOHdLy5csVDAY1efJk1dXVRV/Y0dLSovj4rx7YPfPMM+ro6NBNN93U5ToVFRV64IEHvtn0AIAhy/PvgQ2GtrY2JScn83tg+Fbi98DwbXbO/B4YAADnCgIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJFMBW7dhgzo6OgZ7DKBfRSKR6D+/8847XT4GrBvI79mmAvbAypUaOXKkFi5cONijAP1i8+bNuuyyy6Ifz5o1Szk5Odq8efMgTgX0j4ULFyrlggsG7Pp9Clh1dbVycnKUkJCgvLw8bd++/YzHv/LKK7rkkkuUkJCgK664QrW1tX0aVvr3/60+9thjRAzmbd68WTfddJMOHDjQZf+BAwd00003ETGYtnDhQj322GPq7OwcsHvEOeeclxM2btyo4uJirV27Vnl5eaqqqtIrr7yi3bt3KzU1tdvx27Zt03//93+rsrJSP/nJT/TSSy/p0Ucf1Y4dOzRx4sRe3bOtrU3Jycld9vni43Xok0/k9/u9jA+cEyKRiC7NzdWn//jHaY/JyszUB3/5i3w+XwwnA765jo4OpVxwQbd4hUIhJSUl9d+NnEdTp051paWl0Y8jkYjLzMx0lZWVPR5/8803u+uvv77Lvry8PPeLX/zitPf44osvXCgUim4tLS1OEhsbGxub4e3IkSNek3NGnp5C7OjoUGNjowoKCqL74uPjVVBQoIaGhh7PaWho6HK8JM2cOfO0x0tSZWWlkpOTo9u4ceO8jAkAOAf985//7NfrDfNy8OHDhxWJRJSWltZlf1pamj788MMezwkGgz0eHwwGT3uf8vJylZWVRT8+cuSILrjgArW0tHR7KhFfaWtrU3Z2tvbv39+/D9O/ZVins2ONeod16p1QKKRx48bpvPPO69fregpYrAQCAQUCgW77k5OT+SLphaSkJNapF1ins2ONeod16p34+P594bunq6WkpMjn86m1tbXL/tbWVqWnp/d4Tnp6uqfjAQDoDU8B8/v9ys3NVX19fXRfZ2en6uvrlZ+f3+M5+fn5XY6XpDfffPO0xwMA0Buen0IsKytTSUmJpkyZoqlTp6qqqkrhcFjz5s2TJBUXFysrK0uVlZWSpHvuuUczZszQE088oeuvv141NTX6y1/+onXr1vX6noFAQBUVFT0+rYivsE69wzqdHWvUO6xT7wzUOnn+PTBJeuqpp/TYY48pGAxq8uTJ+u1vf6u8vDxJ0g9/+EPl5OTo+eefjx7/yiuvaOnSpdq3b5++973vafXq1Zo1a1a/fRIAgKGnTwEDAGCwmXovRAAATiFgAACTCBgAwCQCBgAw6ZwJ2GD+iRZLvKzT+vXrNX36dI0dO1Zjx45VQUHBWdf128Dr19IpNTU1iouL0+zZswd2wHOE13U6cuSISktLlZGRoUAgoIsuumhI/HfndZ2qqqp08cUXa8SIEcrOztaCBQv0xRdfxGjawfHOO++osLBQmZmZiouL02uvvXbWc7Zu3aqrrrpKgUBAF154YZdXrvdav741cB/V1NQ4v9/vNmzY4P72t7+5+fPnuzFjxrjW1tYej3/vvfecz+dzq1evdh988IFbunSpGz58uHv//fdjPHlseV2nOXPmuOrqardz5063a9cu9/Of/9wlJye7v//97zGePHa8rtEpe/fudVlZWW769Onupz/9aWyGHURe16m9vd1NmTLFzZo1y7377rtu7969buvWra6pqSnGk8eW13V68cUXXSAQcC+++KLbu3eve+ONN1xGRoZbsGBBjCePrdraWrdkyRK3efNmJ8m9+uqrZzx+z549buTIka6srMx98MEH7sknn3Q+n8/V1dV5uu85EbBY/ImWbwOv6/R1J0+edKNHj3YvvPDCQI046PqyRidPnnRXX321+93vfudKSkqGRMC8rtMzzzzjxo8f7zo6OmI14jnB6zqVlpa6H/3oR132lZWVuWnTpg3onOeS3gRs4cKF7vLLL++yr6ioyM2cOdPTvQb9KcRY/YkW6/qyTl93/PhxnThxot/fEfpc0dc1euihh5SamqrbbrstFmMOur6s0+uvv678/HyVlpYqLS1NEydO1MqVKxWJRGI1dsz1ZZ2uvvpqNTY2Rp9m3LNnj2pra3njhq/pr+/hg/5u9LH6Ey3W9WWdvm7RokXKzMzs9oXzbdGXNXr33Xf17LPPqqmpKQYTnhv6sk579uzRn/70J916662qra1Vc3Oz7r77bp04cUIVFRWxGDvm+rJOc+bM0eHDh3XNNdfIOaeTJ0/qzjvv1P333x+Lkc043ffwtrY2ff755xoxYkSvrjPoj8AQG6tWrVJNTY1effVVJSQkDPY454SjR49q7ty5Wr9+vVJSUgZ7nHNaZ2enUlNTtW7dOuXm5qqoqEhLlizR2rVrB3u0c8rWrVu1cuVKPf3009qxY4c2b96sLVu2aMWKFYM92rfSoD8C40+09E5f1umUxx9/XKtWrdJbb72lK6+8ciDHHFRe1+jjjz/Wvn37VFhYGN3X2dkpSRo2bJh2796tCRMmDOzQg6AvX0sZGRkaPny4fD5fdN+ll16qYDCojo4O+f3+AZ15MPRlnZYtW6a5c+fq9ttvlyRdccUVCofDuuOOO7RkyZJ+/3tYVp3ue3hSUlKvH31J58AjMP5ES+/0ZZ0kafXq1VqxYoXq6uo0ZcqUWIw6aLyu0SWXXKL3339fTU1N0e2GG27Qtddeq6amJmVnZ8dy/Jjpy9fStGnT1NzcHA28JH300UfKyMj4VsZL6ts6HT9+vFukTkXf8bazUf32Pdzb60sGRk1NjQsEAu755593H3zwgbvjjjvcmDFjXDAYdM45N3fuXLd48eLo8e+9954bNmyYe/zxx92uXbtcRUXFkHkZvZd1WrVqlfP7/W7Tpk3uH//4R3Q7evToYH0KA87rGn3dUHkVotd1amlpcaNHj3a//OUv3e7du90f/vAHl5qa6h5++OHB+hRiwus6VVRUuNGjR7vf//73bs+ePe6Pf/yjmzBhgrv55psH61OIiaNHj7qdO3e6nTt3OkluzZo1bufOne6TTz5xzjm3ePFiN3fu3Ojxp15Gf99997ldu3a56upquy+jd865J5980o0bN875/X43depU9+c//zn672bMmOFKSkq6HP/yyy+7iy66yPn9fnf55Ze7LVu2xHjiweFlnS644AInqdtWUVER+8FjyOvX0n8aKgFzzvs6bdu2zeXl5blAIODGjx/vHnnkEXfy5MkYTx17XtbpxIkT7oEHHnATJkxwCQkJLjs72919993uX//6V+wHj6G33367x+81p9ampKTEzZgxo9s5kydPdn6/340fP94999xznu/Ln1MBAJg06D8DAwCgLwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACTCBgAwCQCBgAw6f8BHxr5p0HBFAwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = SinglePhaseFlowWithGeometry()\n",
    "model.prepare_simulation()\n",
    "\n",
    "print(\"The grid:\")\n",
    "pp.plot_grid(model.mdg, plot_2d=True)\n",
    "print(\"One fracture inside:\")\n",
    "model.fracture_network.plot()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, what happenes inside `model.prepare_simulation`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    def prepare_simulation(self) -> None:\n",
      "        \"\"\"Run at the start of simulation. Used for initialization etc.\"\"\"\n",
      "        # Set the material and geometry of the problem. The geometry method must be\n",
      "        # implemented in a ModelGeometry class.\n",
      "        self.set_materials()\n",
      "        self.set_geometry()\n",
      "\n",
      "        # Exporter initialization must be done after grid creation,\n",
      "        # but prior to data initialization.\n",
      "        self.initialize_data_saving()\n",
      "\n",
      "        # Set variables, constitutive relations, discretizations and equations.\n",
      "        # Order of operations is important here.\n",
      "        self.set_equation_system_manager()\n",
      "        self.create_variables()\n",
      "        self.initial_condition()\n",
      "        self.reset_state_from_file()\n",
      "        self.set_equations()\n",
      "\n",
      "        self.set_discretization_parameters()\n",
      "        self.discretize()\n",
      "        self._initialize_linear_solver()\n",
      "        self.set_nonlinear_discretizations()\n",
      "\n",
      "        # Export initial condition\n",
      "        self.save_data_time_step()\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(inspect.getsource(model.prepare_simulation))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `prepare_simulation` method runs a sequence of initialization steps, all of which are important. \n",
    "However, we will focus on two steps for now:\n",
    "\n",
    "- Initializing the independent variables of the simulation in the `create_variables` method. In our case, these are pressure $P$ and interface flux.\n",
    "- Initializing the equations of the model in the `set_equations` method. In our case, these are the mass balance equation and the interface equations.\n",
    "\n",
    "The bookkeeping of variables and equations is done with the `equation_system` attribute of the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "EquationSystem for mixed-dimensional grid with 2 subdomains and 1 interfaces.\n",
       "Variables present on at least one grid or interface:\n",
       "\tpressure, interface_darcy_flux\n",
       "In total 3 equations, with names: \n",
       "\tmass_balance_equation, interface_darcy_flux_equation, well_flux_equation"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.equation_system"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variables creation\n",
    "The variables are created in the method `create_variables`. Inside this method, we register all the model variables in the `equation_system`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    def create_variables(self) -> None:\n",
      "        \"\"\"Assign primary variables to subdomains and interfaces of the\n",
      "        mixed-dimensional grid.\n",
      "\n",
      "        \"\"\"\n",
      "        self.equation_system.create_variables(\n",
      "            self.pressure_variable,\n",
      "            subdomains=self.mdg.subdomains(),\n",
      "            tags={\"si_units\": \"Pa\"},\n",
      "        )\n",
      "        # Note that `interface_darcy_flux_variable` is not multiplied by rho * mu^-1.\n",
      "        # However, after multiplication, whe know that the resulting flux should be a\n",
      "        # mass flux with units  `kg * s^-1`. The units of `interface_darcy_flux` can\n",
      "        # then be inferred by solving the below equation for `int_flux_units`:\n",
      "        # kg * s^-1 = [kg * (m^nd)^-1] * [Pa * s]^-1 * intf_flux_units\n",
      "        self.equation_system.create_variables(\n",
      "            self.interface_darcy_flux_variable,\n",
      "            interfaces=self.mdg.interfaces(codim=1),\n",
      "            tags={\"si_units\": f\"m^{self.nd} * Pa\"},\n",
      "        )\n",
      "        self.equation_system.create_variables(\n",
      "            self.well_flux_variable,\n",
      "            interfaces=self.mdg.interfaces(codim=2),\n",
      "            tags={\"si_units\": f\"m^{self.nd} * Pa\"},\n",
      "        )\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(inspect.getsource(model.create_variables))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Within PorePy, it is a common practice to store the variable names as an attribute of the model.\n",
    "The name of, for instance, the pressure variable can be obtained like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'pressure'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.pressure_variable"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accessing the variable itself can be done by using the variable name like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Mixed-dimensional variable with name pressure, id 28\n",
       "Composed of 2 variables\n",
       "Total size: 20"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = model.equation_system.md_variable(model.pressure_variable)\n",
    "p"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike the `AdArray`, a variable is an abstract AD object that does not store any values.\n",
    "It is used to define equations to be evaluated later. \n",
    "However, the model stores the current state of the values for our variables, and we can access them using the `value` method. It returns a numerical type (a scalar, a Numpy array or a sparse matrix):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_initial = p.value(model.equation_system)\n",
    "p_initial[:10]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have not run the simulation yet, it returns the initial condition which is zero everywhere by default\n",
    "\n",
    "We can change the values manually by calling the `set_variable_values` method of the `equation_system`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.        ,  0.84147098,  0.90929743,  0.14112001, -0.7568025 ,\n",
       "       -0.95892427, -0.2794155 ,  0.6569866 ,  0.98935825,  0.41211849])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "p_new = np.sin(np.arange(p_initial.size))\n",
    "\n",
    "model.equation_system.set_variable_values(\n",
    "    values=p_new,\n",
    "    variables=[p],\n",
    "    iterate_index=0,  #   | For a more advanced reader:\n",
    "    time_step_index=0,  # | We reference method documentation to see what these keyword\n",
    "    additive=False,  #    | arguments do.\n",
    ")\n",
    "\n",
    "# Check that it worked\n",
    "p.value(model.equation_system)[:10]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting equations\n",
    "\n",
    "### A simplified example\n",
    "To set the equations in the model, we use the `set_equation` method inside `prepare_simulation`. This method can be a bit complicated, so we will start with a simplified example.\n",
    "\n",
    "\n",
    "To create an equation, we apply arithmetic operators or functions to variables. \n",
    "For example, we will define a polynomial equation using the pressure variable we saw in the previous section:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_eq = p**2 + p * 5 + 2"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Note for an advanced reader:</b> \n",
    "The variable object incorporates the variables on all the grids in the mixed-dimensional\n",
    "setting. This means, we can operate with one object and simultainously define an \n",
    "equation on all the grids.\n",
    "Here, we will define an equation both on 1D and 2D grids.\n",
    "\n",
    "\n",
    "Since the equation is an AD object, just like the variable, it also has the `value` method.\n",
    "When we run it, the `poly_eq` equation uses the state of `p` to return the AD array and Jacobian matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.        ,  6.91542834,  7.37330894,  2.7255149 , -1.21126246,\n",
       "       -1.87508561,  0.68099553,  5.71656438,  7.92562097,  4.23043407])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poly_eq.value(model.equation_system)[:10]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main purpose of the AD framework is to automatically compute the partial derivatives matrix (Jacobian) of an operator. This can be done with the `value_and_jacobian` method. Note that this method is more computationally expensive than `values`, so use it only when you need the Jacobian. The method returns an instance of `AdArray`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ad array of size 20\n",
       "Jacobian is of size (20, 28) and has 20 elements"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poly_ad = poly_eq.value_and_jacobian(model.equation_system)\n",
    "poly_ad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[5.        , 0.        , 0.        , 0.        , 0.        ],\n",
       "        [0.        , 6.68294197, 0.        , 0.        , 0.        ],\n",
       "        [0.        , 0.        , 6.81859485, 0.        , 0.        ],\n",
       "        [0.        , 0.        , 0.        , 5.28224002, 0.        ],\n",
       "        [0.        , 0.        , 0.        , 0.        , 3.48639501]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poly_ad.jac.todense()[:5, :5]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Note for an advanced reader:</b> \n",
    "In order to make our equation a full-fledged PDE, we must discretize spatial derivatives. \n",
    "This is a more advanced topic, so we refer you to the corresponding tutorials on [flux](./flux_discretizations.ipynb) and [stress](./stress_discretization.ipynb) discretizations."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting equations in PorePy\n",
    "The `set_equations` method is where equations are defined in the model, just like we did with the polynomial equation.\n",
    "For convenience, separate parts of an equation are defined in separate functions. These functions are called separately, so that we assemble an equation from smaller blocks.\n",
    "\n",
    "As an example, let's look at the mass balance equation.\n",
    "This equation is typically defined using the accumulation, flux and source terms. \n",
    "Specific functions for assembling each of these terms can be called separately, before they later are combined to create the full mass balance equation.\n",
    "\n",
    "Below we show how this is done in practice by inspecting the source code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    def mass_balance_equation(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:\n",
      "        \"\"\"Mass balance equation for subdomains.\n",
      "\n",
      "        Parameters:\n",
      "            subdomains: List of subdomains.\n",
      "\n",
      "        Returns:\n",
      "            Operator representing the mass balance equation.\n",
      "\n",
      "        \"\"\"\n",
      "        # Assemble the terms of the mass balance equation.\n",
      "        accumulation = self.fluid_mass(subdomains)\n",
      "        flux = self.fluid_flux(subdomains)\n",
      "        source = self.fluid_source(subdomains)\n",
      "\n",
      "        # Feed the terms to the general balance equation method.\n",
      "        eq = self.balance_equation(subdomains, accumulation, flux, source, dim=1)\n",
      "        eq.set_name(\"mass_balance_equation\")\n",
      "        return eq\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(inspect.getsource(model.mass_balance_equation))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For illustrating how each of these terms are separate from the others, we will access the flux term and evaluate it.\n",
    "The Jacobian we achieve corresponds to that of the flux term only:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ad array of size 49\n",
       "Jacobian is of size (49, 28) and has 72 elements"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subdomains = model.mdg.subdomains()\n",
    "fluid_flux = model.fluid_flux(subdomains)\n",
    "fluid_flux_data = fluid_flux.value_and_jacobian(model.equation_system)\n",
    "fluid_flux_data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Going back to the source code shown above, you can see we use the method `balance_equation` to define \n",
    "a conservation law.\n",
    "Inspecting this method as well shows us that this is where the different parts of the mass balance equation come together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    def balance_equation(\n",
      "        self,\n",
      "        subdomains: list[pp.Grid],\n",
      "        accumulation: pp.ad.Operator,\n",
      "        surface_term: pp.ad.Operator,\n",
      "        source: pp.ad.Operator,\n",
      "        dim: int,\n",
      "    ) -> pp.ad.Operator:\n",
      "        \"\"\"Balance equation that combines an accumulation and a surface term.\n",
      "\n",
      "        The balance equation is given by\n",
      "        .. math::\n",
      "            d_t(accumulation) + div(surface_term) - source = 0.\n",
      "\n",
      "        Parameters:\n",
      "            subdomains: List of subdomains where the balance equation is defined.\n",
      "            accumulation: Operator for the cell-wise accumulation term, integrated over\n",
      "                the cells of the subdomains.\n",
      "            surface_term: Operator for the surface term (e.g. flux, stress), integrated\n",
      "                over the faces of the subdomains.\n",
      "            source: Operator for the source term, integrated over the cells of the\n",
      "                subdomains.\n",
      "            dim: Spatial dimension of the balance equation.\n",
      "\n",
      "        Returns:\n",
      "            Operator for the balance equation.\n",
      "\n",
      "        \"\"\"\n",
      "\n",
      "        dt_operator = pp.ad.time_derivatives.dt\n",
      "        dt = self.ad_time_step\n",
      "        div = pp.ad.Divergence(subdomains, dim=dim)\n",
      "        return dt_operator(accumulation, dt) + div @ surface_term - source\n",
      "\n"
     ]
    }
   ],
   "source": [
    "accumulation = model.fluid_mass(subdomains)\n",
    "flux = model.fluid_flux(subdomains)\n",
    "source = model.fluid_source(subdomains)\n",
    "\n",
    "eq = model.balance_equation(\n",
    "    subdomains=subdomains,\n",
    "    accumulation=accumulation,\n",
    "    surface_term=flux,\n",
    "    source=source,\n",
    "    dim=1,  # Scalar equation\n",
    ")\n",
    "\n",
    "print(inspect.getsource(model.balance_equation))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The complete set of equations is registred in the `equation_system`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    def set_equations(self):\n",
      "        \"\"\"Set the equations for the mass balance problem.\n",
      "\n",
      "        A mass balance equation is set for all subdomains and a Darcy-type flux relation\n",
      "        is set for all interfaces of codimension one.\n",
      "\n",
      "        \"\"\"\n",
      "        subdomains = self.mdg.subdomains()\n",
      "        codim_1_interfaces = self.mdg.interfaces(codim=1)\n",
      "        # TODO: If wells are integrated for nd=2 models, consider refactoring sorting of\n",
      "        # interfaces into method returning either \"normal\" or well interfaces.\n",
      "        codim_2_interfaces = self.mdg.interfaces(codim=2)\n",
      "        sd_eq = self.mass_balance_equation(subdomains)\n",
      "        intf_eq = self.interface_darcy_flux_equation(codim_1_interfaces)\n",
      "        well_eq = self.well_flux_equation(codim_2_interfaces)\n",
      "        self.equation_system.set_equation(sd_eq, subdomains, {\"cells\": 1})\n",
      "        self.equation_system.set_equation(intf_eq, codim_1_interfaces, {\"cells\": 1})\n",
      "        self.equation_system.set_equation(well_eq, codim_2_interfaces, {\"cells\": 1})\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(inspect.getsource(model.set_equations))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During the simulation, the values of the variables will update, and the AD equations will automaticaly produce the Jacobian based on the new values.\n",
    "\n",
    "The full Jacobian and the values of the equations can be accessed through the `equation_system`. This is done inside `pp.run_time_dependent_model`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "jacobian, vals = model.equation_system.assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.08466849, -4.32480821, -3.7754852 ,  0.71956398,  3.81444645,\n",
       "        3.3735598 ,  1.86050941, -4.42321869, -5.82804816, -0.78575445])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vals[:10]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we can see where the values within the full Jacobian is located using the spy function within matplotlib:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fdc4d9562c0>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGiCAYAAACh/hJSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAa5ElEQVR4nO3dfWxT5/nG8StQ4kKbmIWQOBmGBdrCVkqmMcgiWtaKCMgkxNsk+jIJKgSChWpAu1ZMLS/bpGxUQlUrVv4arFKBDqmAijQkCE1Qt8AEBSG0NSJZNoJIQosUO4RiEHl+f0R1f4ZAEmP7Pna+H+mo+PjE5/bjE189OY9vZznnnAAASLEh1gUAAAYnAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABgIm0CaPv27fre976nhx9+WGVlZfrnP/9pXVLKbd68WVlZWTHLpEmTrMtKuuPHj2vevHkqLi5WVlaWDhw4EHO/c04bN25UUVGRhg8froqKCl24cMGm2CTqaxyWLVt21/Exd+5cm2KTpLq6WtOmTVNOTo4KCgq0YMECNTQ0xGxz48YNVVVVadSoUXr00Ue1ePFitbe3G1WcHP0Zh2efffau42HVqlVGFfcuLQLoo48+0vr167Vp0yZ9/vnnKi0t1Zw5c3TlyhXr0lLuySefVGtra3T57LPPrEtKuq6uLpWWlmr79u293r9161a9++672rFjh06ePKlHHnlEc+bM0Y0bN1JcaXL1NQ6SNHfu3JjjY8+ePSmsMPnq6upUVVWlEydO6MiRI7p165Zmz56trq6u6Dbr1q3TJ598on379qmurk6XL1/WokWLDKtOvP6MgyStWLEi5njYunWrUcX34NLA9OnTXVVVVfT27du3XXFxsauurjasKvU2bdrkSktLrcswJcnt378/eru7u9sFAgH39ttvR9d1dHQ4n8/n9uzZY1Bhatw5Ds45t3TpUjd//nyTeqxcuXLFSXJ1dXXOuZ7XftiwYW7fvn3Rbf797387Sa6+vt6qzKS7cxycc+6nP/2p+9WvfmVXVD94/gzo5s2bOn36tCoqKqLrhgwZooqKCtXX1xtWZuPChQsqLi7W+PHj9dJLL+nixYvWJZlqbm5WW1tbzPHh9/tVVlY2KI+P2tpaFRQUaOLEiVq9erWuXr1qXVJShUIhSVJeXp4k6fTp07p161bM8TBp0iSNHTs2o4+HO8fhGx9++KHy8/M1efJkbdiwQdevX7co754esi6gL1999ZVu376twsLCmPWFhYX64osvjKqyUVZWpl27dmnixIlqbW3Vli1b9Mwzz+j8+fPKycmxLs9EW1ubJPV6fHxz32Axd+5cLVq0SCUlJWpqatJvfvMbVVZWqr6+XkOHDrUuL+G6u7u1du1azZgxQ5MnT5bUczxkZ2dr5MiRMdtm8vHQ2zhI0osvvqhx48apuLhY586d0xtvvKGGhgZ9/PHHhtXG8nwA4VuVlZXRf0+ZMkVlZWUaN26c/vrXv2r58uWGlcELnn/++ei/n3rqKU2ZMkUTJkxQbW2tZs2aZVhZclRVVen8+fOD4jro/dxrHFauXBn991NPPaWioiLNmjVLTU1NmjBhQqrL7JXn/wSXn5+voUOH3jWLpb29XYFAwKgqbxg5cqSeeOIJNTY2Wpdi5ptjgOPjbuPHj1d+fn5GHh9r1qzRoUOH9Omnn2rMmDHR9YFAQDdv3lRHR0fM9pl6PNxrHHpTVlYmSZ46HjwfQNnZ2Zo6dapqamqi67q7u1VTU6Py8nLDyuxdu3ZNTU1NKioqsi7FTElJiQKBQMzxEQ6HdfLkyUF/fFy6dElXr17NqOPDOac1a9Zo//79OnbsmEpKSmLunzp1qoYNGxZzPDQ0NOjixYsZdTz0NQ69OXv2rCR563iwngXRH3v37nU+n8/t2rXL/etf/3IrV650I0eOdG1tbdalpdSrr77qamtrXXNzs/v73//uKioqXH5+vrty5Yp1aUnV2dnpzpw5486cOeMkuW3btrkzZ864//3vf8455/7whz+4kSNHuoMHD7pz5865+fPnu5KSEvf1118bV55Y9xuHzs5O99prr7n6+nrX3Nzsjh496n70ox+5xx9/3N24ccO69IRZvXq18/v9rra21rW2tkaX69evR7dZtWqVGzt2rDt27Jg7deqUKy8vd+Xl5YZVJ15f49DY2Oh++9vfulOnTrnm5mZ38OBBN378eDdz5kzjymOlRQA559x7773nxo4d67Kzs9306dPdiRMnrEtKuSVLlriioiKXnZ3tvvvd77olS5a4xsZG67KS7tNPP3WS7lqWLl3qnOuZiv3WW2+5wsJC5/P53KxZs1xDQ4Nt0Ulwv3G4fv26mz17ths9erQbNmyYGzdunFuxYkXG/U9ab89fktu5c2d0m6+//tr98pe/dN/5znfciBEj3MKFC11ra6td0UnQ1zhcvHjRzZw50+Xl5Tmfz+cee+wx9+tf/9qFQiHbwu+Q5ZxzqTvfAgCgh+evAQEAMhMBBAAwQQABAEwQQAAAEwQQAMAEAQQAMEEAAQBMpFUARSIRbd68WZFIxLoUU4xDD8ahB+PQg3HokU7jkFYfRA2Hw/L7/QqFQsrNzbUuxwzj0INx6ME49GAceqTTOKTVGRAAIHMQQAAAE577Qrru7m5dvnxZOTk5ysrKirkvHA7H/HewYhx6MA49GIcejEMPL4yDc06dnZ0qLi7WkCH3Ps/x3DWgS5cuKRgMWpcBAHhALS0t9/2iPM+dAeXk5KRsX6FQKGX78vv9A/6ZVNYHAIkSDocVDAb7fD9PWgBt375db7/9ttra2lRaWqr33ntP06dP7/Pn7vyzWzJ5fYaI1+sDgPvp6/08KZMQPvroI61fv16bNm3S559/rtLSUs2ZM0dXrlxJxu4AAGkoKQG0bds2rVixQi+//LJ+8IMfaMeOHRoxYoT+/Oc/J2N3AIA0lPAAunnzpk6fPq2KiopvdzJkiCoqKlRfX3/X9pFIROFwOGYBAGS+hAfQV199pdu3b6uwsDBmfWFhodra2u7avrq6Wn6/P7owAw4ABgfzD6Ju2LBBoVAourS0tFiXBABIgYTPgsvPz9fQoUPV3t4es769vV2BQOCu7X0+n3w+X6LLAAB4XMLPgLKzszV16lTV1NRE13V3d6umpkbl5eWJ3h0AIE0l5XNA69ev19KlS/XjH/9Y06dP1zvvvKOuri69/PLLydgdACANJSWAlixZoi+//FIbN25UW1ubfvjDH+rw4cN3TUwAAAxenusF9813WaSCx576XeLpCuH15wQg/cTboaav7yQynwUHABicCCAAgAkCCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmEhKN+xE6KuJXSLE22AvHvE0CY3nZ2hgCtjKxN/BgdbX36bSnAEBAEwQQAAAEwQQAMAEAQQAMEEAAQBMEEAAABMEEADABAEEADBBAAEATBBAAAATBBAAwAQBBAAwQQABAEx4ths24pOqDtrx7gvIdPxe9B9nQAAAEwQQAMAEAQQAMEEAAQBMEEAAABMEEADABAEEADBBAAEATBBAAAATBBAAwAQBBAAwQQABAEx4thmp3+8f0PbxNABMZdPAeBt+DlQqxyGe50SjRgDf4AwIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACc82Iw2FQsrNze339jTGTL14xo/XCbDlpd9BzoAAACYIIACAiYQH0ObNm5WVlRWzTJo0KdG7AQCkuaRcA3ryySd19OjRb3fykGcvNQEAjCQlGR566CEFAoFkPDQAIEMk5RrQhQsXVFxcrPHjx+ull17SxYsX77ltJBJROByOWQAAmS/hAVRWVqZdu3bp8OHDev/999Xc3KxnnnlGnZ2dvW5fXV0tv98fXYLBYKJLAgB4UJZL8ocsOjo6NG7cOG3btk3Lly+/6/5IJKJIJBK9HQ6HFQwGM+5zQPHUFw+vf2bG668TkOlS8TsYDofl9/v7fB9P+uyAkSNH6oknnlBjY2Ov9/t8Pvl8vmSXAQDwmKR/DujatWtqampSUVFRsncFAEgjCQ+g1157TXV1dfrvf/+rf/zjH1q4cKGGDh2qF154IdG7AgCksYT/Ce7SpUt64YUXdPXqVY0ePVpPP/20Tpw4odGjRyd6VwCANJbwANq7d2+iH7JfUtUYM959pepCeqomO0ipGwcmLsBKJh57XqqPXnAAABMEEADABAEEADBBAAEATBBAAAATBBAAwAQBBAAwQQABAEwQQAAAEwQQAMAEAQQAMEEAAQBMJP0L6bws3qZ8mdig0Mu83mgWmYvjIbk4AwIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmBjU3bDjlaruzPHsJ5Xde+PtOD1QqRwHOp0DqcMZEADABAEEADBBAAEATBBAAAATBBAAwAQBBAAwQQABAEwQQAAAEwQQAMAEAQQAMEEAAQBMEEAAABM0I00RLzcwxbd4nZAu4m0G7KVjjzMgAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJmhG6mFeb4yZqqaG8TZdjEc8z8nrrxMyo3HnnbxcW39xBgQAMEEAAQBMDDiAjh8/rnnz5qm4uFhZWVk6cOBAzP3OOW3cuFFFRUUaPny4KioqdOHChUTVCwDIEAMOoK6uLpWWlmr79u293r9161a9++672rFjh06ePKlHHnlEc+bM0Y0bNx64WABA5hjwJITKykpVVlb2ep9zTu+8847efPNNzZ8/X5L0wQcfqLCwUAcOHNDzzz//YNUCADJGQq8BNTc3q62tTRUVFdF1fr9fZWVlqq+v7/VnIpGIwuFwzAIAyHwJDaC2tjZJUmFhYcz6wsLC6H13qq6ult/vjy7BYDCRJQEAPMp8FtyGDRsUCoWiS0tLi3VJAIAUSGgABQIBSVJ7e3vM+vb29uh9d/L5fMrNzY1ZAACZL6EBVFJSokAgoJqamui6cDiskydPqry8PJG7AgCkuQHPgrt27ZoaGxujt5ubm3X27Fnl5eVp7NixWrt2rX7/+9/r8ccfV0lJid566y0VFxdrwYIFiawbAJDmBhxAp06d0nPPPRe9vX79eknS0qVLtWvXLr3++uvq6urSypUr1dHRoaefflqHDx/Www8/nLiqAQBpL8t5rKNdOByW3+9XKBTielCKeL1Ro9ebkaYKDUyRLvr7Pm4+Cw4AMDgRQAAAEwQQAMAEAQQAMEEAAQBMEEAAABMEEADABAEEADBBAAEATBBAAAATBBAAwAQBBAAwQQABAEwM+OsYkHni7Zicqu7MqezonKrO214fByAVOAMCAJgggAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABggmakiFs8zTFT1cAUSCeD9feCMyAAgAkCCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmCCAAAAmaEaKlPJ6A9NUNXiM5znFIxMaVt4pExt3er2+ZOEMCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAmakcLzUtXANN59IbV4jTIHZ0AAABMEEADAxIAD6Pjx45o3b56Ki4uVlZWlAwcOxNy/bNkyZWVlxSxz585NVL0AgAwx4ADq6upSaWmptm/ffs9t5s6dq9bW1uiyZ8+eByoSAJB5BjwJobKyUpWVlffdxufzKRAIxF0UACDzJeUaUG1trQoKCjRx4kStXr1aV69evee2kUhE4XA4ZgEAZL6EB9DcuXP1wQcfqKamRn/84x9VV1enyspK3b59u9ftq6ur5ff7o0swGEx0SQAAD8pyDzCpPisrS/v379eCBQvuuc1//vMfTZgwQUePHtWsWbPuuj8SiSgSiURvh8NhBYNBhUIh5ebmxlsaBjmvfw4o3voGis/MwEI4HJbf7+/zfTzp07DHjx+v/Px8NTY29nq/z+dTbm5uzAIAyHxJD6BLly7p6tWrKioqSvauAABpZMCz4K5duxZzNtPc3KyzZ88qLy9PeXl52rJlixYvXqxAIKCmpia9/vrreuyxxzRnzpyEFg4ASG8DDqBTp07pueeei95ev369JGnp0qV6//33de7cOf3lL39RR0eHiouLNXv2bP3ud7+Tz+dLXNUAgLT3QJMQkqG/F68A9C1Vkx0kb0/g8Njb3F0y7Tl5ZhICAAC9IYAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYGPDXMQBAb1LV0Tmen/H6V7Rn4nPqD86AAAAmCCAAgAkCCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmKAZKZDB4m08GW+jy1TsJ1UNTKXU1RePTHhOnAEBAEwQQAAAEwQQAMAEAQQAMEEAAQBMEEAAABMEEADABAEEADBBAAEATBBAAAATBBAAwAQBBAAwQTNSAHeJp/lkqhqYplKqxiFVDUzj3VeyXlvOgAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJigGSlgIFUNK+NtIhnPvlLVUDMTm3163UCfUzgclt/v73M7zoAAACYIIACAiQEFUHV1taZNm6acnBwVFBRowYIFamhoiNnmxo0bqqqq0qhRo/Too49q8eLFam9vT2jRAID0N6AAqqurU1VVlU6cOKEjR47o1q1bmj17trq6uqLbrFu3Tp988on27dunuro6Xb58WYsWLUp44QCA9JblHuCK2ZdffqmCggLV1dVp5syZCoVCGj16tHbv3q2f//znkqQvvvhC3//+91VfX6+f/OQnfT7mNxevQqGQcnNz4y0N8LRMnISQKl6fhID+v48/0DWgUCgkScrLy5MknT59Wrdu3VJFRUV0m0mTJmns2LGqr6/v9TEikYjC4XDMAgDIfHEHUHd3t9auXasZM2Zo8uTJkqS2tjZlZ2dr5MiRMdsWFhaqra2t18eprq6W3++PLsFgMN6SAABpJO4Aqqqq0vnz57V3794HKmDDhg0KhULRpaWl5YEeDwCQHuL6IOqaNWt06NAhHT9+XGPGjImuDwQCunnzpjo6OmLOgtrb2xUIBHp9LJ/PJ5/PF08ZAIA0NqAzIOec1qxZo/379+vYsWMqKSmJuX/q1KkaNmyYampqousaGhp08eJFlZeXJ6ZiAEBGGNAZUFVVlXbv3q2DBw8qJycnel3H7/dr+PDh8vv9Wr58udavX6+8vDzl5ubqlVdeUXl5eb9mwAEABo8BTcO+1/THnTt3atmyZZJ6Poj66quvas+ePYpEIpozZ47+9Kc/3fNPcHdiGjYGA6Zhx49p2N7X3/fxB/ocUDIQQLgTbzipl2ljnolB7GUp+RwQAADxIoAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYiOsbUeFdmdbFWPJ+fZkoVV/9kKrXNt79ePk5ZQLOgAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJgY1M1I42k0KHm72aCXa0Nmy7QGpvHuy+vPyUs4AwIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJgggAIAJAggAYIIAAgCYIIAAACYIIACACQIIAGBiUDcjHawNAJFYXm4+6fWGu5n4O0gD0/7jDAgAYIIAAgCYIIAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJgggAICJQd2MFKmXiU0XvVxfvLVl4uvkZalqYBrvvpKFMyAAgAkCCABgYkABVF1drWnTpiknJ0cFBQVasGCBGhoaYrZ59tlnlZWVFbOsWrUqoUUDANLfgAKorq5OVVVVOnHihI4cOaJbt25p9uzZ6urqitluxYoVam1tjS5bt25NaNEAgPQ3oEkIhw8fjrm9a9cuFRQU6PTp05o5c2Z0/YgRIxQIBBJTIQAgIz3QNaBQKCRJysvLi1n/4YcfKj8/X5MnT9aGDRt0/fr1ez5GJBJROByOWQAAmS/uadjd3d1au3atZsyYocmTJ0fXv/jiixo3bpyKi4t17tw5vfHGG2poaNDHH3/c6+NUV1dry5Yt8ZYBAEhTWS7OSeGrV6/W3/72N3322WcaM2bMPbc7duyYZs2apcbGRk2YMOGu+yORiCKRSPR2OBxWMBhUKBRSbm5uPKXBw/h8SXrgdfI+L38OKBwOy+/39/k+HtcZ0Jo1a3To0CEdP378vuEjSWVlZZJ0zwDy+Xzy+XzxlAEASGMDCiDnnF555RXt379ftbW1Kikp6fNnzp49K0kqKiqKq0AAQGYaUABVVVVp9+7dOnjwoHJyctTW1iZJ8vv9Gj58uJqamrR792797Gc/06hRo3Tu3DmtW7dOM2fO1JQpU5LyBAAA6WlA14Du9TfHnTt3atmyZWppadEvfvELnT9/Xl1dXQoGg1q4cKHefPPNfl/P6e/fDpGeuLaQHnidvG/QXQPqq/BgMKi6urqBPCQAYJCiG7aHZeL/hXq9PvRIVXdmjof4ZUKnc5qRAgBMEEAAABMEEADABAEEADBBAAEATBBAAAATBBAAwAQBBAAwQQABAEwQQAAAEwQQAMAEAQQAMJExzUi91GAvUbxeH/D/0cA0PaTqdeoPzoAAACYIIACACQIIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJgggAIAJAggAYMJzveC+6VMUDoeTvq9U7ANAYvF7mz766jvnuQDq7OyUJAWDwaTvy+/3J30fABKL39v00dnZed/XK8t5rLVsd3e3Ll++rJycnLs6sIbDYQWDQbW0tCg3N9eoQnuMQw/GoQfj0INx6OGFcXDOqbOzU8XFxRoy5N5Xejx3BjRkyBCNGTPmvtvk5uYO6gPsG4xDD8ahB+PQg3HoYT0O/TlTZRICAMAEAQQAMJFWAeTz+bRp0yb5fD7rUkwxDj0Yhx6MQw/GoUc6jYPnJiEAAAaHtDoDAgBkDgIIAGCCAAIAmCCAAAAmCCAAgAkCCABgggACAJgggAAAJv4P6a52jho7T+8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.spy(jacobian.toarray())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What we have explored\n",
    "We have looked at the AD framework and how it is utilized by the PorePy models.\n",
    "\n",
    "Key messages:\n",
    "- The AD framework is capable of computing the Jacobian of an equation automatically. \n",
    "- In PorePy model, the AD objects are prepared in the `prepare_simulation` method.\n",
    "- The variables are defined in the `create_variables` method.\n",
    "- The equations are assembled from the variables in the `set_equations` method.\n",
    "- The process of assembling an equation is modular. You can access separate terms of an equation by calling the model methods.\n",
    "- The equations can be evaluated many times with the new values of the variables:\n",
    "    - with `values` if you need only values.\n",
    "    - with `values_and_jacobian`, and the Jacobian will be computed automatically."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix: AD functions\n",
    "\n",
    "This is the recipe for how to define an AD function. We will use our pressure variable from above and define the `sin` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Mixed-dimensional variable with name pressure, id 28\n",
       "Composed of 2 variables\n",
       "Total size: 20"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AD Operator function 'sin_function'"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sin = pp.ad.Function(func=pp.ad.functions.sin, name=\"sin_function\")\n",
    "sin"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`func` can be anything that takes `AdArray` and returns `AdArray`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Operator with no name formed by Operations.evaluate with 2 children."
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = sin(p)\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ad array of size 20\n",
       "Jacobian is of size (20, 28) and has 20 elements"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = sin(p).value_and_jacobian(model.equation_system)\n",
    "result"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
